#include "Astronomy/RiseTransitSet.hpp"

#include "Astronomy/Celestial/Moon.hpp"
#include "Astronomy/Celestial/Earth.hpp"
#include "Astronomy/CDynamicalTime.hpp"
#include "Astronomy/Precession.hpp"

// 地平线;
#define RadHorizon Transformation::ArcsecToRadians(-50.0 * 60.0)
// 民用;
#define RadCivil Transformation::ArcsecToRadians(-06.0 * 3600.0)
// 航海;
#define RadNautical Transformation::ArcsecToRadians(-12.0 * 3600.0)
// 天文;
#define RadAstronomical Transformation::ArcsecToRadians(-18.0 * 3600.0)

double RiseTransitSet::CalcHourAngle(double h, double w)
{
    double c = (sin(h) - sin(this->Lat) * sin(w)) / cos(this->Lat) / cos(w);
    if (abs(c) > 1)
        return Constants::PI;
    return acos(c);
}

void RiseTransitSet::Mcoord(double jd, bool H0, Twilight &r)
{
    // 低精度月亮赤经纬;
    auto z = CelestialMoon::Coordinate((jd + this->dt) / 36525.0, 40, 30, 8);
    // 转为赤道坐标;
    z.EclipticToEquatorial(this->E);
    // 得到此刻天体时角;
    r.H = Transformation::MapToNegPIToPIRange(CDynamicalTime::pGST(jd, this->dt) + this->Lon - z.X);

    // 升起对应的时角;
    if (H0)
        r.SetHourAngle(Twilight::Type::Horizon, this->CalcHourAngle(0.7275 * Constants::EarthRadius / z.Z - Transformation::ArcsecToRadians( 34.0 * 60.0), z.Y));
}

Twilight RiseTransitSet::Mt(double jd)
{
    this->dt = CDynamicalTime::TDUT(jd);
    this->E  = Precession::EclipticObliquity(jd / 36525.0);

    double n;
    // 查找最靠近当日中午的月上中天,mod2的第1参数为本地时角近似值;
    jd -= modf(0.1726222 + 0.966136808032357 * jd - 0.0366 * this->dt + this->Lon / Constants::PI2, &n);
    const double sv = Constants::PI2 * 0.966;

    Twilight r(jd);
    // 月球坐标;
    this->Mcoord(jd, true, r);
    r.CalcEntry(Twilight::Type::Horizon, sv);

    const auto current = &r.GetEntry(Twilight::Type::Horizon);

    r.Rise             = current->Dawn;
    r.Set              = current->Dusk;

    r.Transit += (0 - r.H) / sv;
    r.LowerTransit += (Constants::PI - r.H) / sv;

    this->Mcoord(r.Rise, true, r);
    r.Rise += Transformation::MapToNegPIToPIRange(-current->HourAngle - r.H) / sv;

    this->Mcoord(r.Set, true, r);
    r.Set += Transformation::MapToNegPIToPIRange(+current->HourAngle - r.H) / sv;

    this->Mcoord(r.Transit, false, r);
    r.Transit += Transformation::MapToNegPIToPIRange(0 - r.H) / sv;

    this->Mcoord(r.LowerTransit, false, r);
    r.LowerTransit += Transformation::MapToNegPIToPIRange(Constants::PI - r.H) / sv;

    return r;
}

void RiseTransitSet::Scoord(double jd, int xm, Twilight &r)
{
    // 太阳坐标(修正了光行差);
    Coordinate3D z{CelestialEarth::EclipticLongitude((jd + this->dt) / 36525.0, 5) + Constants::PI - Transformation::ArcsecToRadians(20.5), 0, 1};
    // 转为赤道坐标;
    z.EclipticToEquatorial(this->E);
    r.H = Transformation::MapToNegPIToPIRange(CDynamicalTime::pGST(jd, this->dt) + this->Lon - z.X);

    // 地平以下50分;
    if ((xm == 10) || (xm == 1))
        r.SetHourAngle(Twilight::Type::Horizon, this->CalcHourAngle(RadHorizon, z.Y));
    // 地平以下6度;
    if ((xm == 10) || (xm == 2))
        r.SetHourAngle(Twilight::Type::Civil, this->CalcHourAngle(RadCivil, z.Y));
    // 地平以下12度;
    if ((xm == 10) || (xm == 3))
        r.SetHourAngle(Twilight::Type::Nautical, this->CalcHourAngle(RadNautical, z.Y));
    // 地平以下18度;
    if ((xm == 10) || (xm == 4))
        r.SetHourAngle(Twilight::Type::Astronomical, this->CalcHourAngle(RadAstronomical, z.Y));
}

Twilight RiseTransitSet::St(double jd)
{
    this->dt = CDynamicalTime::TDUT(jd);
    this->E  = Precession::EclipticObliquity(jd / 36525.0);

    double n;
    // 查找最靠近当日中午的月上中天,mod2的第1参数为本地时角近似值;
    double y = modf(jd + this->Lon / Constants::PI2, &n);
    jd -= y;

    const double sv = Constants::PI2;

    Twilight r(jd);
    // 太阳坐标;
    this->Scoord(jd, 10, r);
    //  地平线;
    r.CalcEntry(Twilight::Type::Horizon, sv);
    // 民用;
    r.CalcEntry(Twilight::Type::Civil, sv);
    // 航海;
    r.CalcEntry(Twilight::Type::Nautical, sv);
    // 天文;
    r.CalcEntry(Twilight::Type::Astronomical, sv);

    // 中天;
    r.Transit += (0 - r.H) / sv;
    // 下中天;
    r.LowerTransit += (Constants::PI - r.H) / sv;

    this->StEntry(Twilight::Type::Horizon,1, r);
    this->StEntry(Twilight::Type::Civil,2, r);
    this->StEntry(Twilight::Type::Nautical,3, r);
    this->StEntry(Twilight::Type::Astronomical,4, r);

    this->Scoord(r.Transit, 0, r);
    r.Transit += Transformation::MapToNegPIToPIRange(0 - r.H) / sv;

    this->Scoord(r.LowerTransit, 0, r);
    r.LowerTransit += Transformation::MapToNegPIToPIRange(Constants::PI - r.H) / sv;

    auto rr =  r.GetEntry(Twilight::Type::Horizon);
    r.Rise = rr.Dawn;
    r.Set  = rr.Dusk;

    return r;
}